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We present an analysis of well-posedness of constrained evolution of 3+1 formulations of GR. In 
this analysis we explicitly take into account the energy and momentum constraints as well as possible 
algebraic constraints on the evolution of high-frequency perturbations of solutions of Einstein's equa- 
tions. In this respect, our approach is principally different from standard analyses of well-posedness 
of free evolution in general relativity. Our study reveals the existence of subsets of the linearized 
Einstein's equations that control the well-posedness of constrained evolution. It is demonstrated 
that the well-posedness of ADM, BSSN and other 3+1 formulations derived from ADM by adding 
combinations of constraints to the right-hand-side of ADM and/or by linear transformation of the 
dynamical ADM variables depends entirely on the properties of the gauge. For certain classes of 
gauges we formulate conditions for well-posedness of constrained evolution. This provides a new 
basis for constructing stable numerical integration schemes for a classical Arnowitt-Deser-Misner 
(ADM) and many other 3+1 formulations of general relativity. 

PACS numbers: 04.25.Dm, 04.70.Bw 



I. INTRODUCTION 

An outstanding problem of numerical general relativity 
(GR) is achieving long-term stable numerical integration 
of Einstein's equations. Until very recently, problems 
such as the general case of colliding black holes (BH) 
could not be solved due to instability of numerical in- 
tegration [l], @, H, 0, [E @1- During the last year, sev- 
eral groups have succeeded in simulating certain binary 
black hole (BBH) spacetimes. Collisions of non-rotating 
black holes were simulated using BSSN f|f with a specific 
choice of "1 + log" slicing and modified gamma- freezing 
shift conditions in Q and A three-dimensional colli- 
sion of two black holes originated from the collapse of a 
scalar field was simulated using a generalized harmonic 
decomposition of the GR field equations in 0]. In addi- 
tion to using a special gauge, the authors of Q and @ 
enforced some of the BSSN constraints and used some of 
the constraints to control the constraint violating modes. 
The author of [{|, in addition to special coordinates he 
used the constraints to damp constraint violating modes. 

In spite of this remarkable success the general problem 
of long-term and stable integration of the Einstein equa- 
tions remains unsolved. There is no general method to 
choose an appropriate formulation and appropriate coor- 
dinates to guarantee well-posedness and stability of nu- 
merical integration for a general strong field case. In 
particular, one does not know how the approaches used 
in 0, H, Q will behave in other strong vacuum field as- 
trophysical cases, or cases where matter or even matter 
and magnetic fields are present. Perhaps the final word to 
this problem would be the development of schemes which 
implement fully constrained evolution of well-posed for- 
mulations. The purpose of this paper is to formulate 
a general approach to study the well-posedness of con- 



strained evolution of 3+1 formulations of GR. 

Generally, a 3+1 formulation is comprised of the evolu- 
tionary part and the constraints. The standard approach 
to solving a 3+1 system consists of integration of the evo- 
lutionary part in time (free evolution) starting with con- 
strained initial conditions. If the constraints are satisfied 
initially, they should automatically be satisfied through- 
out the evolution due to the mathematical properties of 
Einstein's equations. 

Well-posedness of free evolution of 3+1 formulations 
has been analyzed in 0, fl5| . For example, a classi- 
cal ADM 3+1 formulation [l| is usually ill-posed. Ill— 
posedness of free evolution precludes stable numerical in- 
tegration. We note, however, that well-posedness can not 
guarantee global in time existence of solutions, but only 
local existence. However, it is a necessary condition for 
stable integration. 

There has been a number of attempts to overcome the 
instability of numerical integration. Using a harmonic 
gauge makes the ADM 3+1 system well-posed [lOfl. How- 
ever, this gauge is not convenient for many physical prob- 
lems. Introduction of a conformal factor and the trace of 
extrinsic curvature as additional unknown variables into 
the system allows to increase the duration of stable inte- 
gration The evolutionary part of a 3+1 system can 
also be modified by adding a combination of constraints 
to its right-hand side. Choosing a special gauge (densi- 
tized lapse and zero shift) and addition of certain com- 
binations of constraints to the right hand sides (RHS) of 
the ADM evolution equations makes the modified system 
well-posed However, all these modifications did not 
lead to a complete elimination of instabilities. Numeri- 
cal experiments show that in general three-dimensional 
problems of GR the constraint equations are eventually 
violated even for a well-posed free evolution, and this ter- 
minates computations. Little progress has been made to- 
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wards a theoretical understandin g of this behavior. Pos- 
sible explanations are given in [6lJl5j. 

Recently, attempts have been made to improve the be- 
havior of modified 3+1 systems by enforcing the con- 
straints after every integration time step of a free evo- 
lution 0, || . According to [|| this procedure allows in- 
tegration of an isolated spherical black hole space-time 
for extended periods of time. We also mention that for 
certain cases perturbative approaches provide an alter- 
native to straightforward numerical integration, e.g., for 
forming initial conditions for BH collisions [l2| . Yet, the 
general problem of long-term stable integration remains 
open. 

In a high-frequency perturbation analysis of a free evo- 
lution it is possible to separate perturbations on three 
parts: (1) perturbations of space-time itself, (2) pertur- 
bations of a coordinate system, and (3) perturbations 
describing deviations from constraints. If the behavior of 
space-time at a given point does not depend on future, 
then we must associate ill-posedness with coordinate and 
constraint- violating modes of perturbations. To achieve 
stable numerical integration, we must (A) use a gauge 
that does not lead to ill-posedness, and (B) eliminate 
or suppress ill-posedness caused by constraint violating 
modes. 

A general theory of gauge stability (problem A) has 
been formulated, and well-posedness of gauges has been 
analyzed in [l3l |. It was demonstrated that coordinate 
perturbations in the evolution of the metric can be sepa- 
rated from the other two types of perturbations and the 
study of gauge stability can be reduced to a study of a 
general quasi-linear system of eight coupled partial dif- 
ferential equations for perturbations of lapse a, shift Pi, 
i = 1, ...,3, and perturbations of space-time coordinates 
x a , a = 0, 3. Conditions for well-posedness have been 
formulated in [13j for several types of gauges. We will 
repeatedly return to this subject in subsequent sections. 

Constraint-violating modes of perturbations (problem 
B) are not fully understood at present. Recently, at- 
tempts have been made to stabilize numerical integration 
by enforcing the constraint equations after every integra- 
tion time step of a hyperbolic free evolution [J H[ . Such 
enforcement is not a unique procedure. Several possi- 
bilities are discussed in According to 0], constraint 
enforcement improves integration of an isolated spheri- 
cal black hole space-time. Analysis of well-posedness of 
constraint enforcing procedure for a hyperbolic 3+1 for- 
mulation, densitized lapse, zero shift, and flat Minkowski 
space-time is given in [111 ]. 

An alternative to enforcement of constraints after a 
free evolution time step may be the construction of nu- 
merical schemes for constrained evolution in which grow- 
ing constraint- violating modes are explicitly removed. In 
order to achieve this goal we must understand the nature 
of evolution of perturbations which satisfy constraints. 

In this paper we present a general analysis of 
constraint-satisfying perturbations and address the issue 
of well-posedness of constrained evolution of 3+1 for- 



mulations of GR. We explicitly take into consideration 
the energy and momentum constraints on the evolution 
of high-frequency perturbations of solutions of Einstein's 
equations. In this respect, our analysis is principally dif- 
ferent from standard analyses of well-posedness of a free 
evolution in general relativity. Our study reveals the ex- 
istence of subsets of the linearized version of Einstein's 
equations that control well-posedness of constrained evo- 
lution. We demonstrate that the well-posedness of ADM, 
BSSN and other 3+1 formulations derived from ADM 
by adding combinations of constraints to the right-hand- 
side (RHS) of ADM and/or by linear transformation of 
the dynamical ADM variables depends entirely on the 
properties of the gauge. For certain classes of gauges 
we formulate conditions of well-posedness. The existence 
of these subsets provides a basis for constructing stable 
numerical integration schemes that incorporate the con- 
straints directly. 

The paper is organized as follows. We begin with 
the ADM 3+1 formulation and gauge classification (Sec- 
tion II) . In Section III we give a general theory of well- 
posedness of constrained evolution of ADM. In Section IV 
we extend our theory to other 3+1 formulations including 
the Kidder-Scheel-Teukolsky (KST) and the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN). Discussion and con- 
clusions are given in Section V. 

II. ADM 3+1 FORMULATION AND GAUGES 

A general form of the ADM 3+1 formulation consists 
of the evolutionary part 

^ = ~2aK lJ + Vi/3j + Vj/3i, (1) 

^ =a ^R tJ + KKij - 2 7 mn K„„J^„) - ViVj-a 
+ {Vi(3 m )K mj + (VjDK^ + p m V m K l3 , 

(2) 

and the energy and momentum constraints which we will 
call the kinematic constraints, 

U : <®R + K 2 - K mn K mn = 0, (3) 

M: \7 m K m i -\7 i K = 0, z - 1,2,3, (4) 

where K = j mn K mn , 7^ and ify are the three- 
dimensional metric of a space-like hypersurface and the 
extrinsic curvature respectively, a is the lapse func- 
tion, P l is the shift vector (these are gauge functions), 
and ( 3 'Rij is the three-dimensional Ricci tensor, ^'R = 
j^^Rij. We must add a specification of gauge (lapse 
and shift) in order to close the system {TJ, 

In this paper we use a general gauge specification sim- 
ilar to that introduced in [l3j with one modification. In- 
stead of working with the dual shift vector (3^ here we 
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work with the shift vector j3 k = 7 fej '/3j. A general gauge 
can be specified as 



F a I x b ,a,0 l 



da d 2 a 



W dx b dx c '"'' dx bl "' lljl dx b ' ' 
a, b,c = 0, ...,3, i, _?' = !, 2,3, 



07* 



= 0, 



(5) 



where it is implied by the ellipsis that one can add higher 
order derivatives of both the lapse and the shift vec- 



tor, e.g. 



' dx a dx b 



or the dynamical variables, e.g., 



d 2 K zj 
dx a dx b ' 



Q x aQ x b , and so on. Following [13j, we distinguish three 
types of gauges. 

1. Fixed gauges for which both the lapse and shift are 
functions of coordinates t = x° and x % , i = 1, 3 only, 

a = a{t,x l ), (3 k = /?*(*, X*), i, fe = 1,2,3. (6) 

Geodesic slicing a = 1, /3 = is a specific case of a fixed 
gauge. 

2. Algebraic (local) gauges for which both the lapse and 
shift can be expressed as algebraic functions of coordi- 
nates and local values of 7y and its derivatives, 



.7ij 



' dx b 



(3 k = (3 h 



>7*j 



(7) 

3. Differential (non-local) gauges which are defined by 
a set of partial differential equations and which cannot 
be reduced to an algebraic form. Algebraic gauges are 
a subset of differential gauges. We call the differential 
gauges non-local, because the gauge variables are not 
completely defined by the local values of the dynamical 
variables. Differential gauges are governed by differential 
equations and hence the lapse and shift may also depend 
on boundary conditions. 

For fixed and algebraic gauges, the total number of par- 
tial differential equations of the ADM formulation does 
not increase compared to ([I])- dU). For differential gauges, 
a complete ADM formulation will consist of (jTJ) - (J3J) plus 
differential equations describing the gauge. 



III. ANALYSIS OF WELL-POSEDNESS 

We begin with a brief description of our general ap- 
proach to the analysis of well-posedness of constrained 
sets of partial differential equations (PDEs). Let 



d t u = M e (u)d e u + M°(u), €=1,2,3 



(8) 



be a set of n first order quasi linear partial differential 
equations, where u is the column vector of the n unknown 
variables, M l are nxn matrices and M° is a n-component 
column vector. 

The concept of mathematical well-posedness is often 
related to strong hyperbolicity. For system ([5]) we present 
the following theorem from [2 3} without proof. 



a. Theorem 1. The Cauchy problem for a first- 
order system of quasi-linear PDEs ([5]) is well posed if 
and only if the following two conditions hold: 

1. For all unit one forms x>i, all eigenvalues of the char- 
acteristic matrix, M = ViM 1 are purely real. 

2. There is a constant K, and for each Vi, there is a 
transformation S(vi) with 



IsCOI + ls-^oi <k, 



(9) 



such that the transformed matrix 
5(?;i)M(wi)S , ~ 1 (wi) is diagonal. 

b. Definition 1. A first-order system of quasi- 
linear PDEs ijHJ) is called strongly hyperbolic if all condi- 
tions of the theorem above are met. 

c. Definition 2. A first-order system of quasi- 
linear PDEs © is called weakly hyperbolic if it satisfies 
only the first condition of Theorem 1 and does not sat- 
isfy the second condition. Weakly hyperbolic systems are 
ill-posed. 

For a constrained system ©, the dynamical variables 
satisfy a set of m < n quasi linear constraint equations 
of the form: 



C e (u)d e u + C°(u) = 0, €=1,2,3 



(10) 



where C e are m x n matrices and C° is an m-component 
column vector . We call a constrained surface a collection 
of all solutions of ([8]) which satisfy the constraints. If 
evolution starts on the constrained surface it will remain 
on this surface. 

In order to study well-posedness one drops the zeroth- 
order terms M°, freezes the coefficients A 1 of ((SJ) and 
studies the characteristic matrix of the system for all 
frozen in problems. This is equivalent to considering 1) 
high frequency and 2) small amplitude planar pertur- 
bations on the dynamical variables along a line locally 
specified by a unit vector v and parameterized by A, so 
that for an arbitrary function u(x x , x 2 , x 3 ) 



du 
dx l 



du 
l 9A' 



(11) 



where Vi is the dual vector to v l , given by w,; = jijV 3 and 
7y is the 3-metric of a spacelike hypersurface embed- 
ded in the manifold carrying the background solution u 
about which we perturb. The two aforementioned prop- 
erties are the properties of all perturbations considered in 
this paper and those will be implied whenever the word 
"perturbation" is used. For perturbations du on u, com- 
bination of ([8]) and (1 1 1 [) yields 

8 t 5u = M%^ = M^. (12) 

CM OA 

For perturbations of u to remain on the constraint sur- 
face, they must satisfy the linearized constraint equa- 
tions. After linearizing equations (|10[) we obtain 



dSu 



C 



dSu 



0. 



(13) 
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where M and C are the principal matrices of equations 
dH]) and (flT)]) respectively. Equations p^|) is a set of m 
equations for the spatial derivatives of the n unknown 
variables, which in general can be solved for m of the 
n spatial derivatives of variables u. Substitution of (| 13() 
into (fT2|) . leads to a set of q = n — m linear partial dif- 
ferential equations for q of the initial n variables. This is 
schematically given by 



at -Mw) ax > 



(14) 



where A q is a (q x q) matrix. We will refer to (JT3J as 
the minimal set. The solution of (p~4|) completely deter- 
mines the solution of the entire linearized system (fT2f . 
Therefore, the well-posedness of the minimal set deter- 
mines the well-posedness of the entire constrained sys- 
tem. Theorem 1 and Definitions 1,2 apply to (fT4|) . 

We should note here that the description above is 
for first-order systems of PDEs. When one deals with 
second-order systems, then it is straightforward to ob- 
tain the equivalent first-order system of PDEs, by sim- 
ply defining the first-order derivatives as new variables. 
Courant and Hilbert [2l| show that if one derives the 
first-order form in the fashion described above, then the 
totality of solutions of the two systems coincide, for given 
Cauchy data. In addition to that in [22| it is shown that 
the hyperbolic properties of the second-order system and 
its equivalent first-order counterpart are the same. Al- 
though a reduction to a first-order system is not neces- 
sary, it makes it easier to analyze well-posedness, since 
in second-order systems one has to carefully study the 
behavior of the first order terms. 



A. Linearized equations of ADM 

We want to study well-posedness of an ADM 3+1 for- 
mulation (evolutionary part fl},© plus constraints ([3]), 
((3]) ). For the analysis of well-posedness it is convenient 
to rewrite equations (jT]) - (J^J) in first-order form. We 
introduce new variables 



Di4: 



ij;k 



dx k 



(15) 



and drop all low-order terms, which do not contribute 
to the principal part of the equations. In terms of these 
variables, the ADM equations linearized with respect to 
a certain unperturbed solution (not necessarily a flat 
Minkowski spacetime) 



7-tj'i 

can be written as 



Ki. 



at 



= a ( - -d.Ld.j5a 
at \ a 



pi 



/3 k dSK i: 
a dx k 



(16) 



(17) 



(18) 



+ T k J d k 5a + 2K k(l d 3) 5f3 k , 



BSD 



ij;k 



at 



d8K l3 t FdSDij-k 



+ 2d k [ le(l d 3) Sp"]+D l3 . e 



2K, 



1 1 dx l ) 13 dx k 

dS[3 e 



dx k ' 



(19) 

where 5"fij, SK^j, SD^j.^, 5a, and 5(3 m are perturbations 
of (fT6|) . Rjj is the principal part of the Ricci tensor 



( 



05D. 



ik;j 



dSDjk-i 05 Djrt, 



ij;k 



2' \ dx 

1 



dx s 



dx s 



1 mn d5D mn;i 



(20) 



and 



dxJ 



r k 3 ^-l kn (D in;j + D 



(21) 



To close the system equations (|T7| - (|23l) must be sup- 
plemented with the linearized version of gauge equations 

We notice that the evolution equations for the 3-metric 
(|17p are decoupled from the evolution equations for Kij , 
Dij-k and the linearized gauge equations. It is the subsys- 
tem (fT8|) . ([T9|) and the linearized gauge equations which 
determine the well-posedness of the entire system. The 
solution of ([17]) is completely determined by the solution 
of the above subsystem. 

The linearized constraint equations in new variables 
are 



and 



n : r 1 - r m Rl m = o, 



dx s 



dx s 



(22) 



(23) 



where 5f is the Kronecker symbol. The introduction of 
extra variables imposes new linear constraints on the sys- 
tem, which, for perturbations of (|16p . can be written as 
follows: 



85D 



ij;k 



05D 



lj;m 



dx r 



dx k 



(24) 



where (|24|) is derived by use of (jl5|) and the fact that 
partial derivatives commute. 



Analysis of well-posedness of constrained 
evolution and well-posed subsets 



We consider planar perturbations 

Sjij, 5Kij, 51), rJ . . 5a, 5.1' 



(25) 



moving in an arbitrary direction locally specified by a 
unit vector v l , v % vi = 1. Substitution of (|25[) into (fl~7|) 
- (|19|) gives a set of thirty linear PDEs for perturba- 
tions (f2l>l) . Substitution of into the eighteen in num- 
ber PDEs gives twelve independent linear PDEs for 



5 



those perturbations. This means that along the direc- 
tion which we are probing only six components of the 
eighteen SDij-k are independent. The linearized energy 
and momentum constraints (|22[) and (|23[) give four addi- 
tional linear PDEs for perturbations (|25| . None of these 
equations contain time derivatives of (125j) . because the 
constraints do not involve time derivatives of the pertur- 
bations. 

Here we note that in the context of the lst-order formu- 
lation equations (fT5|) are constraints on the initial data. 
They are not involved in the analysis of well-posedness 
but they serve to guarantee the coincidence of the solu- 
tions of the 1st order and second order ADM equations. 
Since the evolution of the perturbations of the 3-metric 
is decoupled from and determined by the evolution of the 
perturbations of the subset variables Kij, Dij.^ and the 
linearized gauge equations, the study of well-posedness 
of a constrained evolution reduces to analyzing the sub- 
set of variables K^, Dij-k, and all remaining constraints 
([2"2]). ([25]). (pM]) . There are twenty four evolution equa- 
tions and sixteen constraints involved that leads to eight 
degrees of freedom. 

By elimination of 16 spatial derivatives of the dynami- 
cal variables from the subset of twenty four equations, we 
obtain eight linearly independent equations, which form 
the minimal set. This can be schematically written as 



da 8 - da 8 
— =A 8 (u,v t ) — , 



(26) 



where A 8 is an 8 x 8 matrix which depends on the direc- 
tion of planar perturbations and the background solution, 
a 8 are the independent perturbation amplitudes. 

To study the well-posedness of the Cauchy problem of 
system (|2l)]) we first find the eigenvalues ojk of the prin- 
cipal matrix of the system for every direction vi . If non- 
zero imaginary parts are present in some of those eigen- 
values, this indicates that part of the system forms an 
elliptic subset and the Cauchy problem is ill-posed. If all 
LUk are real, then we investigate whether A 8 is diagonal- 
izablc for every direction Vk (uniformly diagonalizable), 
and whether the second condition of Theorem 1 is 
satisfied (transformation S is uniformly bounded). 

A convenient way to investigate properties of A 8 is to 
use a reduction of A 8 to a Jordan canonical form [171 ] , 



A, 



sjs- 



so that 



-( 

dt \ 



= J 



(27) 



(28) 



where S is a non-singular matrix of a similarity transfor- 
mation (|2"7]) and J is a block-diagonal Jordan canonical 
matrix 



(29) 



Jl 


. 


. 





h ■ 


. 





. 


■ Jn 



consisting of canonical Jordan blocks . Each canonical 
block Jfc has the form 



Jk 








. 


. 


" 


1 




. 


. 








1 


Wfc . 


. 











. 1 


OJk. 



(30) 



In general, the number of canonical blocks N may vary 
from one to eight and the size of each square block can 
vary from one to eight as well. The diagonal elements 
of J are eigenvalues of A 8 . If for any and all possible 
directions Vi all Jordan blocks have size one (i.e. all off- 
diagonal elements of J are zero) , then A 8 has a complete 
set of eigenvectors and hence A 8 is uniformly diagonaliz- 
able. 

If some of off-diagonal elements of J are non-zero at 
least in one direction Vk, the set of eigenvectors in this 
direction is not complete and hence (|26p cannot admit 
a well-posed Cauchy problem. The system (|26p is then 
weakly hyperbolic. 

In general, uniform diagonalizability of A 8 may not 
be enough to guarantee the existence of a uniformly 
bounded transformation S required by Theorem 1. We 
show in Appendix A, however, that if the eigenvalues 
ujk of A 8 are analytic functions and the elements of A 8 
are ratios of analytic functions, then uniform diagonaliz- 
ability guarantees the existence of a uniformly bounded 
transformation matrix S. For the general relativity cases 
studied in the paper, the elements of the matrices of all 
minimal sets arc always ratios of polynomials in Vk ■ Fur- 
thermore, those matrices have real eigenvalues which are 
analytic functions of Vk- Thus, for those systems stud- 
ied here uniform diagonalizability is sufficient in order to 
satisfy the conditions for well-posedness. 

The Jordan decomposition also provides informa- 
tion about combinations of independent variables that 
evolve according to the corresponding eigenfrequencies 
LUk- These combinations can be determined by applying 
a similarity transformation U~ x to the vector of original 
variables a 8 (|2"8"|) . In the case of week hyperbolicity or 
elliptic behavior, this allows one to find a combination 
of variables whose time evolution is responsible for this 
behavior. 

The above discussion is strictly valid for ADM and 
any 3+1 formulation of general relativity derived from 
ADM by a linear transformation of variables and addi- 
tion of combinations of constraints on the RHS of the 
ADM equations, provided that they are coupled with 
fixed or algebraic gauges. If a gauge is differential (speci- 
fied by a general ©), then the number of linearized ADM 
equations may be greater than thirty and the final sys- 
tem of independent perturbations (|26p will contain more 
than eight components. Analogous consideration must 
then be applied to this larger reduced system. A de- 
tailed discussion of well-posedness of ADM 3+1 coupled 
with differential gauges is given in section IIII Dl Next 
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section discusses well-posedness of ADM with fixed and 
algebraic gauges. 



hyperbolicity) of the constrained ADM 3+1 formulation 
with algebraic gauges (f3"Tj) and (|3"2")l therefore are 



C. Well-posed subsets for ADM with fixed and 
algebraic gauges 

For fixed and algebraic gauges, the study of (j2T>)) can 
be carried out analytically but resulting expressions for 
an arbitrary direction Vk are extremely complicated. We 
present here only results for algebraic gauges of a simple 
form 

a = a{t,x k , lik ), l? = f?(t,x k ) (31) 



and 



a = a(t,x k , 7ik ), 0i = I3.it, x k ). (32) 



For gauge (|31|) the perturbation frequencies obtained 
from (l26l) are 



da2 
V 9 7y 



(33) 



and for gauge (|3"2")) they are 



,..,6=±^ 



A/3 1 , 



^7,8 



The difference in formulas (|3"3")) and (|3~i|) arises because 
(|3Tj) fixes the shift vector whereas (j3"2")l fixes its dual 
counterpart. Metric- independent contravariant shift cor- 
responds to metric-dependent covariant shift and vice- 
versa. The results for gauge (|3"2")) coincide exactly with 
those of our analysis of this gauge in [l3| . 

Six eigen-frequencies ujx,...,6 in (|3"3"| and (|3^|) are 
real. Eigenfrequencies ujr,8 in ([3"3"|) and (|34[) are real, if 

(l^j) ^7 > and (f^J + ftft') v i v i > °> respectively. 
The eigenvectors of A s in (|26p are uniformly linearly inde- 
pendent, for this particular choice of gauges. For the case 
of fixed gauges, — 0, the eigenvalues are still real, 

but ^8 has only seven linearly independent eigenvectors. 
That is, all fixed gauges lead to ill-posed constrained evo- 
lution. The general conditions for well-posedness (strong 



da 2 \ 

- — ViVj > 



(35) 



and 



[da 2 



+ P(3 3 ViVj > 



V»* "77-';- . < 36) 

for an arbitrary Vi , respectively. Again, (|36[) is the result 
obtained in 



The analysis of matrix As in (j2l)|) for algebraic gauges 
shows that the minimal set (|26l) separates into four inde- 
pendent subsets each consisting of two equations. Among 
these subsets there are three well-posed subsets corre- 
sponding to pairs of eigenvalues u>x,...,6' These three sub- 
sets describe propagation of gravitational waves and a 
gauge wave. We call these wave subsets. The fourth 
subset, which corresponds to eigenvalues w 7j8 , describes 
another two gauge modes. The solution of the fourth 
subset depends on solutions of the strongly hyperbolic 
wave subsets and the gauge choice. Those four subsets 
completely describe the behavior of Einstein's equations 
(evolutionary part + constraints) in the high frequency 
limit. It is therefore evident that the posedness of Ein- 
stein's equations depends entirely upon the properties of 
the gauge. 

As a simple illustration, we explicitly present the form 
of these subsets for a general metric 7^ , extrinsic curva- 
ture K%j, gauge a = a('jij,x n ,t) and (3 k = f] k (x n ,t), and 
propagation of perturbations along x 1 coordinate direc- 
tion, Vk — (i>i, 0, 0) and 7 11 i>it>i = 1. The eigenvalues for 
this case are: {/3 l vi — a, /3 1- i>i — a, /3 1- i>i + a, (3 l v\ + a, 
fi-vi, fi-vx, frvx + a^/2A ll v 1 v 1 , (Pvx - a^/2A 11 v 1 v 1 } 
and the minimal set is explicitly given by the following 
four subsets 



d6K 



2:S 



dt 

dSD 23 ;l 



dt 



BSK- 



:s:i 



II 



dt 



dSD 



33;1 



dt 



-avx 



-avx 



-av\ 



-av\ 



1 11 dSD 23 . 1 1 dSK 23 



'2 7 



dX a OX 
dSK 23 , P 1 d5D 23 . 1 



OX 



a 



OX 



(37) 



- 2 



l^d SDss-x P l dSK 33 
dX a dX 
8SK 33 l3 x d8D 33 , 



2 7 



OX 



ox 



(38) 



7 



d&D 12 . x j 2 (7 1 v v 3 + 7 1 W 2 - zyv v 3 ) a^ 23 2 (y w 3 - 7 1 W 3 ) a^ 33 

— a^i < 



dt 



in 



7 11( 7 12 7 12 _ 7 11 7 22) 

/3 1 dSD 12 , 



dX 



7 11( 7 12 7 12 _ 7 11 7 22) 



OX 



dSD 



13;1 



dt 



-avi 



2 



a dX 
12 



(39) 



cW^23 



9A 



+ 7 1 



OX 



+ 



(i 1 0SD 13 . 1 
a dX 



dSK n 
dt 



_ A ndSDn ]1 _ 2A 
dX 



ox 



ox 



- ( §7 d Ai + Y A + 2A 2d + 



<9A 



IV : 



1 22 .33 , 1 33 , ,33 , ,22 . 33 \ gggggg /3 1 dSK n 
2 7 a + 2 7 + + a a dx + dx 



(40) 



dSDna 
dt 



-avi 



OX 



dX 



where 



A 11 = 



j33 



<91n< 

07ij 



d 23 = -2 



7 12 7 13 _ 7 11 7 23 

^yl2ryl2 ,yll,-y22 



7 13 7 13 



7 " 7 33 



(41) 



7 12 7 12 



7 11 7 22- 



As one can easily see these equations are valid only when 
7 12 7 12 _ 7 n 7 22 ^ 0. This has been our assumption to 
solve the momentum and Hamiltonian constraints for the 
derivatives of the dynamical variables which we wanted to 
eliminate. Although it may seem that this is not a general 
result we point out that the kinematic constraints can 
always be used to eliminate 4 of the dynamical variables 
provided that certain conditions are held true. If the 
condition above is not satisfied then there will be another 
set of variables that we will be able to eliminate and thus 
obtain the minimally coupled set of partial differential 
equations. In this respect we have not lost generality. 

Subsets I, II, and III describe wave propagation and 
are well-posed. The first two propagate with the shift 
plus the fundamental speed (the lapse function) and the 
third one with the shift velocity. Subset IV will be well 



posed, if the lapse satisfies 



> 0, 



din a „, - n , da 

-5 VtVT, > U => a 

0711 0711 

which is a particular case of the general condition ([55)1 . 
If -tr— = (fixed gauge) , the fourth subset takes the form 



On, 



85 K n 
dt 



IV: 



85D U; 
dt 



avi 



2 ' 



avi 



2 ' 



7 



23 



dSD 



23:1 



dx 



1 

2 7 



33 



dSD 33 g | ^ dSK^ 
dX a dX 

dSK n {P_d5D n .y 
dX a dX 



(42) 



This subset is weakly hyperbolic and ill-posed. This can 
be most easily seen if we consider a simplest case with 
SD23-A — SD 33 -i = and /3 1 = 0. Then the solution 
of (03]) will be SK n = SK n (X,0) and SD n;1 (X,t) = 
5£>ii ; i(A,Q) + (^^(A.O)) t. The linear growth of 
5Dii-i depends on initial conditions and may be arbi- 
trarily fast. Since in the high frequency limit we can 
treat the gauge functions as constants, the physical ac- 
celeration can be neglected. Therefore, the linear growth 
of <5-Dii ; i physically describes the deformation of a syn- 
chronous reference frame with time, which in a general 
non-linear case when the perturbations are not small, 
leads to the formation of caustics. 

Constrained evolution of perturbations of all other 
variables is completely determined by the solution of sub- 
sets I - IV. As an example we present the evolution of 
8K\ 3 and 5Dz2;i when, for simplicity, the shift vector is 
zero : 



dSK 13 
dt 

dSD 2 2;l 

dt 



-av± 



1 



12 



dSD 2 



13 



dSD. 



33:1 



= — 2avi 



dx 

(7 U 7 23 



12„,13 



7 7 



dx 

) dSK 



23 



( 7 H 7 



33 



( 7 12 7 12 
7 13 7 13) 



7 11 7 22) 



dx 



d5K, 



33 



( 7 12 7 12 



7 H 7 22 ) 



dX 



(43) 

The amplitudes of 5K\ 3 and 8D22A will not grow. Equa- 
tions for 8K22 and 8K12 are analogous but more compli- 
cated. All other perturbations satisfy equations 



dSD. 



y';2 



dSD 



jj;3 



0. 



(44) 



dt dt 

We found that the behavior described above is similar 
to that of a general case of algebraic gauges and any 
arbitrary direction of propagation of perturbations. 



8 



Examples of algebraic gauges are the widely used gauge 
a = C(x 1 )^ 1 ^ 2 often referred to as the "harmonic" gauge 
0, the "1+log" gauge a = 1 + log(7) and the densi- 
tized lapse gauge a = C{x l )^ (J Q, all depending on the 
determinant of the three-metric, 7 = det(jij). For these 
gauges, condition (|35|) can be written as 



<31n< 

91k 



91n a dj 



<91nc 

dj 



■JJ V ViVj 



<91n < 
din 7 



> 0. 

(45) 

It can be readily seen that both "harmonic" and "1+log" 
gauges satisfy this condition and lead to a well posed 
constrained evolution. The densitized lapse will provide 
a well-posed constrained evolution only if 



a > 0. 



(46) 



D. Well-posedness of ADM with differential gauges 

Similar approach to posedness of 3+1 formulations can 
be carried out for more complex gauges involving non- 
zero shift and general elliptic, parabolic, or hyperbolic 
differential gauges. As examples of elliptic and parabolic 
differential gauges we will consider the maximal slicing 
condition and its parabolic extension. Although these 
two gauges are believed to prevent coordinate singular- 
ities, here we demonstrate that both of them produce 
weakly hyperbolic minimal sets and thus produce ill- 
posed constrained evolution. 

The maximal slicing condition [l6[ is K = 0, where K 
is the trace of the extrinsic curvature. In the following 
analysis we need the evolution equations for K and the 
determinant 7 of the 3-metric in vacuum. Those can be 
derived by taking the traces of |T]) and ((2|) and they are 



d t In 7 



1/2 



-olK + Vi/3* 



(47) 



d t K = -7 ij 'ViVja + aK i3 K i] + (3 i V l K (48) 



If the K = condition is imposed then (|48[) results in 
the following elliptic differential equation for the lapse 
function 

7 ij 'V i V j Q - aK i:j K ij = 0. (49) 

If the Hamiltonian constraint ([3]) is satisfied (|49]) yields 

7 li V i V J a - aR = 0. (50) 

In the limit of high frequency perturbations the Ricci 
scalar vanishes on the surface of constraints and equation 
(15*01 reduces to 



^drfjSa = 0, (51) 
which, if written along a given direction v l , yields 
d 2 Sa d 2 Sa 



OX 2 d\ 2 



0. 



(52) 



Then the principal term of WiWjSa — ViVj^g- vanishes 
and therefore the ADM + Maximal Slicing equations 
have the same properties as ADM + Fixed Gauges, which 
means that the Cauchy problem for ADM+Maximal Slic- 
ing is ill posed. We should keep in mind that result ((52|) 
is not only valid for a constrained evolution, but also for 
an unconstrained one, since it can be derived from (|49[) . 
too. Equation (|52|) results because of the high frequency 
perturbations we are considering here. 

Here we ought to resolve the apparent contradiction 
between our analysis and the fact that maximal slicing 
prevents the formation of coordinate singularities [lg]. 
This can be seen if (1471) is written as 



K = £ ft ln( 7 



l/2\ 



(53) 



where £^ denotes the Lie derivative along the unit nor- 
mal vector n to the spacelike hypersurfaces with 3 metric 
jij . If we set K = then from (|53p the local volume ele- 
ment 7 1 / 2 is proper-time independent and cannot shrink 
to zero. This means that a coordinate singularity cannot 
be formed. 

However, the well-posedness properties of algebraic 
K = slicing condition are different from those of (f50|) 
with initial conditions K = 0. If maximal slicing is im- 
posed by using K — at all times by eliminating one 



of the components, for example Kn 



7 'J-<5 1, I 5 1J 7 11 



-Kij 

this reduces the number of dynamical degrees of freedom 
in the minimal set l|2"6"l) from eight to seven [27j • In this 
case the perturbations of the trace of the extrinsic cur- 
vature are identically zero, SK — 0, and the minimal set 
of seven equations is well posed. On the other hand, if 
we use (f5U)) and initial conditions K(t = 0) = alone, 
the minimal set of eight equations is ill-posed because 
SK are not necessarily zero. If we write equation l|48p in 
the limit of high frequency perturbations, keeping (|5"2"|) 
in mind, we obtain 



d t SK = (3 l di5K. 



(54) 



Although this equation is well-posed and hence the per- 
turbations of K do not grow, we notice that the K = 
condition may now be violated by perturbations of K. 
It is this violation which is the root to the ill-posedness 
associated with (f5"0")) . 

Let us illustrate this with the following example. Con- 
sider planar high frequency perturbations along x 1 about 
Minkowski spacetime, which means that the unperturbed 
lapse is a = 1 and the unperturbed shift is (5 l = 0. In 
this case, the minimal set for the differential maximal 
slicing condition contains ([42]) , which for this case can be 
written as 



IV 



OSKn 
dt 
dSDn-i 
dt 



0, 



-2v x 



dSK n 



(55) 



This subset (|55p is weakly hyperbolic and hence ill-posed. 
In addition, if we impose the algebraic maximal slicing 
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condition at all times, then in the high frequency limit 

(56) 



d8K 1± = 8SK 22 dSK 33 



<9A v <9A <9A 

This equation eliminates SKu from the minimal set. The 
linearized momentum constraints for this particular case 
also read 

dSK 2 2 , d8K 33 



dX 



dX 



= 0. 



(57) 



Combining ^ and (J37J) we obtain ^|"> = and (55]) 
reduces to one equation 



rv : -?pi = o. 



(58) 



Equation (|58[) is well posed, and this eliminates the pos- 
sibility of formation of coordinate singularities. 

The parabolic extension of maximal slicing [19| has the 
following form 



da _ 1 

dt ~ e 



-■' l I),l),<\ - KijK v a - cK) , (59) 



where e is a positive constant. Since the lower order 
term —K^K^a — cK does not belong to the principal 
part, application of high frequency perturbations along 
any arbitrary direction yields 



dSa 1 d 2 Sa 
~dt = l~d>?' 



(60) 



The linearized evolution equations for variables belonging 
to the minimal set for this particular gauge are given by 
(see (PJ| and (fT5]0: 



d8K K 
dt 



1 d 2 Sa , \ dSa ,„,. 

-ViVi ^rr^ + RU ) + FiiVk — , (61) 



ij;k 



Of 



dX 2 



OA 



c9A 



2K ll v k 



dSa 
~dX' 



(62) 



In the high frequency limit we are considering here, the 
evolution equation for the lapse (|6H)) is completely decou- 
pled from (|6Tj) and ([62| whatsoever. But, the linearized 
"parabolic maximal slicing" is a diffusion equation, which 
is well-posed and it dictates that the amplitude of the 
perturbation of the lapse function in this case diffuses 
out and hence it does not grow. This means that as time 
passes the entire system will asymptotically resemble the 
case of fixed lapse and zero shift, which is a special case 
of a fixed gauge. Therefore the constrained evolution of 
ADM with "parabolic maximal slicing" is ill-posed. In 
appendix [B] we present a more rigorous proof of this fact. 
These results for maximal slicing and its parabolic exten- 
sion agree with those obtained in [l3| for the same slicing 
conditions. 

Following Bona et al. [14] hyperbolic gauges can be 
given in general by the following conditions: 



dta = — a 2 Q, 
d t p l =-aQ\ 



(63) 



where Q, Q % will be given by either algebraic or differen- 
tial equations relating them with other variables of the 
system and will be chosen accordingly in order to obtain 
hyperbolic equations for the lapse and/or the shift. For 
such slicings, one has to modify the ADM equations by 
defining new variables in order to obtain the 1st order 
form. Therefore we define Ai = c^lna, Bf = djf3 l and 
with these definitions we compute: didjU — u(AiAj + 
diAj) and didj(3 l — diB/. Thus, the linearized principal 
part of the ADM equations can be written as 



dt 

d5K tj 
dt 
dSDjj-k 
dt 



-2^8 Bj) , 



=a ( R] 



2 dSK l3 
dx k 



(3 k dSKjj _ dSAj 
a dx k dx l 
(3 e d8D ij;k 



(64) 



a 
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dSB/ 
dx k 



dx e 

dSB^i 



dx k 



We have introduced twelve new variables and therefore 
we need twelve additional evolution equations to describe 
them. One could expect that the number of differential 
equations of the minimal set would be 8+12=20 equa- 
tions. However, we must impose eight new linear con- 
straints arising due to the introduction of new variables 
B k and Ai, 



and 



dB % k _ dBj k 
dxi dx % 



dAj _ dA 3 
dxi dx % 



(65) 



(66) 



For perturbations moving along direction Vi and assum- 
ing vi we get 



and 



dBj k _ Vi dBi k 
dX vi dX 

dAj = VjdAt 
dX vi dX 



(67) 



(68) 



Equations (|67|) and (|68|) mean that there are 3 linearly 
independent B k and 1 independent Ai. Therefore the 
minimal set will in principle consist of only twelve evolu- 
tion equations for our dynamical variables. 

As an example we will consider the Bona-Masso family 
of slicing conditions [l8| which we write in terms of the 
new variables as 



d t In a = ?Ai - af(a)(K - K a ), 



(69) 



where /(a) is a strictly positive function of the lapse, K is 
the trace of the extrinsic curvature and K a = K(t = t Q ). 
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Then the principal part of the evolution equations of the 
perturbations of the new variables is 



dtSAi = (3 k d k 8Ai - a/(a) 7 



hi 



dSK. 



kl 



dx % 



(70) 



However, because of constraints ([68]) only one of the 
above equations will be a part of the minimal set. Our 
analysis of the minimal set, carried out for this gauge, 
shows that the constrained evolution is well posed. 

As an illustration of the latter, we present the lin- 
earized constrained evolution equations for the zero shift 
vector case (also known as the K-driver condition). The 
minimal set then consists of only nine partial differential 
equations (three equations get eliminated due to j3 l = 0) 
for planar perturbations of the dynamical variables Kn, 



K 23 , K 33 , D U]1 ,Di2;i, D 13 -x, D 23 :i, Aj3;i, and A x along 
the a; 1 -axis, which we group into the following subsets 



asK- 



2:t 



I : 



at 

d5D 231 



dt 



dSK : 



33 



II 



dt 
dSD 33 . A 



dt 



1 n^£> 23; i 
-a«i7 



2avi 



dSK. 



23 



OX 



1 11 d5D 33 , 1 
-awi7 



dX 



= — 2av\- 



dSK : 



33 



OX 



(71) 



(72) 



dSDi 



III 



dt 
dSD 13 , 1 



dt 



; — 2avi 
2avi 



( 7 12 7 12 7 13 + 7 11 7 13 7 22 _ 27 11 7 12 7 23) gg^ ^12^3^3 



7 11 7 12 7 33^) 



dSK, 



33 



12 dSK 23 
7 — +1 



7 11( 7 12 7 12 _ 7 11 7 22) 

13 dSK 33 " 



OX 



7 11( 7 12 7 12 _ 7 11 7 22) 



dX 



dX 



(73) 



IV 



dSK u 
dt 



dSAx 
dt 



1 



22^23 



•7 



23 



dSD 



23; 1 



22^33 



„33' 



dX 

dSD 33 , x dAi 



dX 



OX 



- af(a)vi 



n dSK n 



+ 2 



7 12 ( ' 7 13 7 22 



dX 

7 12 7 23 ) d5K. 



23 



7 12 7 12 



7 11 7 22 Q X 
7 13 7 13 7 22 _ 7 12 7 12 7 33 Qfi^ 

dX~ 



7 11 7 22 



and an equation for the evolution of SDn-i 



dSD u 
dt 



-2v\a. 



dSKu 
dX ' 



(74) 



(75) 



Subsets I and II describe gravitational waves propagat- 
ing with the fundamental speed. The eigenvalues which 
correspond to these two subsets are {—a,— a, a, a}. 
These are well-posed and completely decoupled from the 
rest of the system. The perturbations which correspond 
to subset III are completely determined by the solution 
of the first two subsets and they do not grow. Two zero 
eigenvalues correspond to the third subset and therefore 
it corresponds to static modes. The fourth subset is cou- 
pled to the first two, but it is not completely defined 
by the solution of I and II. Subset IV is also well posed 
and it d escribes a gauge wave propagating with speed 
awfja). Finally equation (|75"|) . which describes a static 
gauge mode, is completely determined by the solution of 
the fourth subset and the perturbation for SDn-i does 



not grow. It is clear therefore that this system of equa- 
tions has a well-posed Cauchy problem. We found exactly 
the same behavior in the most general case of perturba- 
tions of arbitrary direction. The Jordan matrix is always 
diagonal and hence the system is well-posed. If the shift 
vector is fixed, but non vanishing, its presence does not 
affect the well-posedness of the constrained system, since 
the Jordan matrix is still diagonal and the 9 non zero 
real eigenvalues are 

f3*v l: P l v % ^ l v ll (i l v l ± a, fat ± a, L3 l v t ± a^/J{a) (76) 

One can demonstrate, as in the previous section, that 
the solution of the entire linearized system can be re- 
trieved once the minimal set is solved and that its posed- 
ness is completely dependent on the posedness of the min- 
imal set. Thus, we conclude that the Bona-Masso family 
of slicing conditions gives rise to a well-posed constrained 
ADM evolution. 



IV. ANALYSIS OF EXTENDED 3+1 
FORMULATIONS 



The analysis of the standard ADM 3+1 formulation 
presented in the previous section can be applied to other 
3+1 formulations of GR. In this section we will explic- 
itly study two of those re-formulations of GR, namely 
the Kidder-Scheel-Teukolsky (KST) (and all KST- 
like formulations) and the Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) @. 
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A. The Kidder-Scheel-Teukolsky formulation of 
3+1 GR 

Kidder et al. 3f suggested a new formulation of Ein- 
stein's equations with a strongly (or even symmetric) hy- 
perbolic set of evolution equations. They obtained this 
formulation by adding terms proportional to the con- 
straint equations to the RHS of the ADM evolution equa- 
tions. This does not change the physics the equations 
describe but changes the character of partial differential 
equations which describe the free evolution. The mod- 
ified set they suggested is (using their notation for this 
section only) 

dtKn =(. . .) + jN 9ij C + CNg ah C a(l])bl 
dtdkij =(• . •) + r)Ng k (iCj) + xNg i3 C k , 

where (...) stands for the RHS of the ADM evolution 
equations, N is the lapse function, g a b is the 3-mctric, 
Kij the extrinsic curvature, dkij are the same variables 
as our Dij-k , C and Ci are the hamiltonian and momentum 
constraints, respectively and 

Ckiij = d[f.dqij = 0. (78) 

Finally, {7, £, 77, x} are arbitrary constants. 

The system was closed with the densitized lapse 

Q = ln(Ng-°), (79) 

where Q is a function independent of the dynamical 
fields, g is the determinant of the three-metric and a is 
the densitization parameter, and was found that a > is 
essential for obtaining a well-posed set of evolution equa- 
tions. This is exactly what we obtain by our analysis 
(|46|) . without adding constraints to the ADM equations, 
but explicitly imposing them. 

If we apply our constrained perturbation analysis to 
the KST formulation we will cancel the added constraints 
on the RHS of their formulation. As result the KST 
formulation has exactly the same analysis as the ADM 
formulation. 

According to [3j any transformation of dynamical vari- 
ables does not change the hyperbolic classification of a 
set of PDEs, if this transformation satisfies the following 
conditions 

1. The transformation is linear in all dynamical vari- 
ables except possible the metric 

2. The transformation is invertible 

3. Time and space derivatives of the metric can be 
written as a sum of only the non-principal terms. 

Their redefinition of variables and the introduction of 
"kinematical" ones is a transformation which satisfies 
the aforementioned criteria and thus it does not affect 



the hyperbolic properties of the set of evolution equa- 
tions. Hence, the constrained perturbation analysis of 
the KST formulation and all KST-like formulations (i.e. 
all those formulations which are derived by addition of 
constraints to the RHS of ADM and perhaps a transfor- 
mation of variables with the aforementioned properties) 
is equivalent to that of the ADM equations. 

This argument may also be used to conclude that any 
3+1 system directly obtained from ADM using the above 
transformation will be equivalent to ADM. 



B. The Baumgarte-Shapiro-Shibata-Nakamura 
formulation of 3+1 GR 

The BSSN formulation was initially introduced by 
Nakamura et al. Q , then modified by Shibata and Naka- 
mura 0], and it was later reintroduced slightly modified 
by T. Baumgarte and S. Shapiro 0. Before we proceed 
with our analysis let us review the BSSN formulation 
first. In what follows, we will use the notation introduced 
by Baumgarte and Shapiro (2j. 



1. Basic variables and equations 

The fundamental dynamical variables of BSSN are 
(ip, ^ij ,K ,Aij ,Y l ) instead of (7^,^), where 





=(l/12)log(det 7 ^), 


lij 


=e" 4v 7 i j, 


K 




Aij 


=e- i *(K ij - (l/3) 7ij K), 


r 





The "connection" symbols r* fe are the Christoffell sym- 
bols associated with the conformal three- metric 7^ . In 
the BSSN formulation, the Ricci curvature tensor is cal- 
culated as 

nBSSN 1 n 

ij — ij ^iji 

Rf 3 = - IDiPw - 2% J D k D k p 

+ l(D t <p)(D 3 ip) - 4n fij (D k <p)(D k ip), 
Ra - - (l/2)7 ,fc fta fc 7y + 7fc (i 9 i) f fe 

+ r r (iJ -) k + 27 m r i( -r i)fem + 7 m r im r fcy , 

(81) 

where Di is the covariant derivative associated with 7^ . 
The evolution equations for these dynamical variables 
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can be written as 



dt<p = - (l/6)aK + (1/6)^(9^) + (d^), 
8t% = - 2aA ij + jik(djp k ) + 7jfe(di/3 fc ) 

-(2/3)%(d k f3 k )+p k (d k j lj ), 
d t K = - D l Dia + aAijAv + (l/3)aK 2 + p(diK), 

o 7 -itpfn n \TF I -4w / U BSSN\TF 

d t Aij=-e ^(DiDja) +e ^ a\R l3 ) 
+ aKA i:j - 2aA ik A k j + (dif3 k )A kj 
+ (d 3 p k )A kt - (2/3)(d k p k )A l3 + p k (d k A lj ), 
d t f l = - 2(d ja )A» + 2a(f i jk A k i - {2/3)f j (djK) 
+ 6A»(d^)) - d 3 ([3 k (d k fi) - t 3 (d k p l ) 
-l kl (d k p) + (2/3)^{d k [3 k )), 

(82) 

where all trace free (TF) two index quantities TJy are 
given by 



T- 



(TF) 



T, 



1 



y -- 7y T, T = ~? kl T kl . (83) 



The constraint equations are 



n BSSN = R BSSN + R 2_ R R ij = q 



M 



BSSN 



= M 



ADM 



= 0, 



g l = r - 7 Jfc r; fc = o, 



A 

s 



Aijfi = o, 
7-1 = 0. 



(84) 
(85) 
(86) 
(87) 
(88) 



f^BSSN anc j j^bssn are Hamiltonian and momen- 
tum constraints (the kinematic constraints) respectively, 
while the latter three are algebraic constraints due to the 
requirements of BSSN formulation. 



2. Linearized equations of BSSN 

The formulation, as given above, is first order in time 
and second order in spatial derivatives. The second or- 
der derivatives occur in the evolution equations for the 
conformal traceless extrinsic curvature. In order to apply 
our analysis we need the first order form. Therefore, we 
define new variables 



d m ip = 4> m and 9 m 7y = D ij;r , 



(89) 



where ";" does not imply a covariant derivative, but it 
separates indices of different nature. In terms of those 
variables the BSSN equations linearized with respect to 
a certain background spacetime solution (not necessarily 
a Minkowski spacetime) 



A, 



K, 



D 



ij;kt 



(90) 



d t 5% =j ik (d j 6l3 k )+j jk {d i 6/3 k ) ~ (2/3)j ij (d k 6(3 k ), 
d t SK = - e-^f k d t d k 5a + P^SK), 
d t 8A l3 = - e-^{didi6a) TF + e -^a(5R* SSN ) TF 
+ (di5p k )A kj + {d 3 5[3 k )A kl 
- (2/3)(d k S(3 k )A l3 + (3 k (d k 5A l3 ), 
d t 6f l = - 2{dj8a)A ij - {A/3)a^ {d 3 5K) 
-f3 k d j (5D ij . >k )+j kj (d j d k 5f3 i ) 
+ j ki (d 3 d k SP) ~ (2/3)fi(d 3 d k 5f3 k ), 

(91) 



d t 5D i3 . k = - 2a - J + f3 l 



dx k ' H dx k ' ^ ii dx^>dx k 
d 2 Sf3 e 2 _ d 2 8(3 l _ - d5a 



8 2 Si3 e 



lit 



dx l dx k 3^ 1 ^ dx e dx k 11 dx k 



2Ai 



Wk 



^(d e 6f3 e )D l3[k , 
1 dSK 1 ■ dS<pi d 2 8j3 i 

I (3 — "I r—— 

6 dx k 6 dx k dx l dx k 

6 dx k 6 dx k ' 
where 

6<p, S%, SA l3 , 6K, Sf\ SD l3 ., k , 54>i (93) 



(92) 



are small amplitude and high frequency perturbations of 
QB9 and 



6R 



BSSN 



= - 2d 3 5(j> l - 2j ij j ek d e Scf >k - ^ K diSD iJik 



( 94 ) 

In the context of the first order formulation, since the 
evolution of the conformal 3-metric is decoupled from the 
evolution of the rest of the system (just like in ADM), 
then the derivative of constraint (|88|) provides a con- 
straint for the Dij- k variables, which has to be taken into 
consideration and is 

f j D ij; k = (95) 
Then, the linearized constraint equations in new variables 



§R BSSN = ^SRBSSN = Q) 



.udSAki 2dSK 
7 --z—r = 0, 



dx l 



3 dx l 



(96) 



(97) 
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d5r ^_ 4m dSD km]j 



dx s 



~ik ~ jm 
= ry ry J 



dx s 



T 



dx k 



(98) 



(99) 



The equations for the constrained evolution of BSSN 
are extremely complicated. For a general background 
solution, we were able to carry out an analytic analysis 
of BSSN with fixed gauges only. Similarly to ADM, we 
found that the constrained evolution of BSSN with fixed 
gauges is ill-posed. 
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dx s 



0. 



(100) 



The introduction of additional variables implies the in- 
troduction of new linear constraint equations which for 
perturbations of ([90|) can be written as follows 



dmSDij-k = d k 5D l 



and d„ 



and 



dmSjij = and d m 5tp = 0. 



(101) 



(102) 



Finally equations ([9"T j) -(f92" l) have to be supplemented 
with the linearized gauge equations (|5j). 



Analysis of well posedness of BSSN with fixed and 
algebraic gauges 



Constraints (|101|) dictate that there is only one inde- 
pendent (j>k and 6 independent L>ij- k . Equation (| 102[) tells 
us (just like in the case of the ADM formulation) that 
the evolution of Sip and 8%j is decoupled from the evo- 
lution of the perturbations of the remaining dynamical 
variables. For BSSN the Hamiltonian constraint can al- 
ways be solved for the derivative of that 4> k involved in it, 
thus eliminating <p k from the minimal set. The momen- 
tum constraints can be solved for the spatial derivatives 
of two Ay's and the spatial derivative of K. Constraint 
(|87|) can be used for the elimination of one of the compo- 



nents of Aij. Finally, constraints (|9"8"|) completely elimi- 
nate the T l variables and (|100[) can eliminate one more 
of the Dij-k variables. This means that fully imposing 
the linearized constraint equations results in a set of lin- 
ear PDEs for three A y - and five Dij-k, the well-posedness 
of which will determine the well-posedness of the entire 
linearized system (|9Tj) - (j9"2"|) . 

In this section we mainly analyze gauges for which the 
lapse function is dependent on the coordinates and the 
dynamical variables and the shift vector is function of 
only the coordinates. With this in mind, the linearized 
principal part of the set of variables which will be part 
of the BSSN minimal set can be written as follows 



d t SD 



ij;k 



e-^aSR?/ SN + [3 k d k 5A lJ1 
2ad k 5A ij + ^dtSb^.u. 



(103) 



To obtain the minimal set one has to fully impose con- 
straints JHHD- PH0|) on equations (fTU3"|) . 



For algebraic gauges, the simplest case possible is per- 
turbations about a flat space 7^ — Sij . If we define 
A y = f^- 21 and A = our analysis shows that the 

minimal set has eigenvalues {(3 l Vi, P" l v il (3 z Vi ± a,f3 l Vi ± 
a,(3 l Vi ± ^r}, where 



D = A + 12A lj v 



4A km S km > 



(104) 



is the necessary condition for all eigenvalues to be real. 
Since the Jordan matrix for this case is diagonal the 
BSSN constrained evolution will be well-posed if (|104|) 
is satisfied. As in the ADM case one could show that 
solving the minimal set of BSSN is adequate to obtain 
the solution of the entire linearized system. 

As an illustration we present a set of constrained evo- 
lution equations for perturbations of variables 



Mn, SA22, SA 23 , 6D U; i, 

<5-Di2;l, 5£>i 3; i, (5£>22;1, ^23:1 



(105) 



propagating along the x 1 direction and perturbed about 
flat space. 



dSA 



2:1 



dt 



1 dSD 
2 a —A 



23;1 



x d5A 



2:! 



dt 



2a 



8SA 



2.'! 



■0 



OA 

1 dSD23a 



(106) 



dSD 12 -i ol dSD 12 -i 

=P 



II 



dt 
d5D 13 . tl 
dt 



--0 



ox 

id SD 13; i 
dX 



(107) 
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dSA u 
dt 



(A" _ A 33 + ri) 

d6D 12 ;l 



A.dSD na 



8 ' d\ 



dX 

(A 22 - A 33 ) 

:dSD 2 3;l 



ox 



22 x 2,3\ ^ SD 22;l 

SA 



III 



9M2 

at 



f 2A 23 — — 
OA 

(A 11 - A 33 + 
<9A 

. (A 22 - A 33 - 

:dSD 2 3;l 



P 1 85 A 



li 



a dX 



A 3 dSDu-i 
¥ ~ 4^ dA 



,<9<SL> 



+ 2A1 dX 

3 dSD 22;1 

2 ] 



13;1 



2A 



23 ' 



9A 



<9A 

/3 1 06 A 



22 



a OX 



dt 
d5D 22 -i 
dt 



= -2a 



2a — + (3 ^x— 

d5A 22 al dSD 22; i 



t)X 



ox 



(108) 



Subset I corresponds to a gravitational wave propagat- 
ing with the shift+ fundamental speed along x 1 . It is de- 
coupled from the other subsets and it is well-posed. Sub- 
set II describes propagation of two waves, which travel 
with the shift vector speed. This subset just like the 
first one is decoupled from the rest of the system. It is 
a wave subset and therefore it is well-posed. Subset III 
is coupled to the first two subsets and hence its solution 
depends on the solutions of I and II. At a first glance 
it may seem as a contradiction that there are 3 subsets 
when for ADM we had 4. However, subset III consists of 
two independent. One of them describes a gravitational 
wave travelling with the shift+fundamcntal speed and 

the other is a gauge wave travelling with speed (3 1 ± 

if D > 0, where for this case 



D = A-4(A n +A 22 + A 33 ) + 12A n , (109) 



which is a special case of (|104[) . 

The constrained evolution of BSSN with fixed gauges, 
that is, A = A u = has the same subsets, but III is now 
decoupled from I and II: 



85 An = gid6Au 



III : 



dt 
d5A 22 
dt 

ggjW 

dt 
ddD 22 -i 
dt 



2a 



= -2a 



dX ' 
a d5Dn ; i 
4 dX 
dSAu 



dX 
dSA 22 
dX 



ad8D 22a X 85A 22 
2 dX +l ~d\~' 



dX ' 

1 d SD 22 - 1 
dX ' 



(110) 



The four eigenfrequencies of this system are {/3 1 , 1 , (3 1 ± 
a}, which are all real. However, the Jordan matrix is 
not diagonal, which implies that the principal matrix of 
(|110| does not have a complete set of eigenvectors, hence 
the constrained evolution is weakly hyperbolic, that is, 
ill-posed. To assign a physical meaning to weak hyper- 
bolicity it is instructive to consider the case of vanishing 
shift. Then one can easily see that subset III breaks into 
two subsets. The ill-posed subset is: 



dSA u 
dt 



= 0, 



dSD Ui i „ dSA u 
= —2a- 



dt 



ox 



(in) 



and its solution is 6 An — SAn(X,t — 0), and 

( 2*|iifA.nV\ t. The linear 



5£>n ; i(A,t) = ^ii ; i(A,0) + (^u(A,0)J t. 

growth of (5Z?ii ; i depends on initial conditions and may 
be arbitrarily fast. As in the ADM formulation, phys- 
ically the linear growth of SDn-i describes the inertial 
deformation of a synchronous reference frame with time 
which, in general non-linear case when perturbations are 
not small leads, to the formation of caustics. 

A special class of algebraic gauges is that with the 
lapse function depending on the determinant of the 3- 
metric, that is, a — 0(7), which in the BSSN formu- 
lation obtains the form a = a(e 12v ), and therefore the 
lapse does not depend on the conformal 3-metric, but 
on the conformal factor only. Thus the strong hyper- 
bolicity condition (|104[) reduces to D = A > 0. For 
"harmonic" slicing a — c(x l )e &v and "1 + log" slicing 
a = 1 + 12(^3 already discussed in the previous section 
it is easy to show that they produce a well-posed con- 
strained BSSN formulation. Similarly, for a densitized 
lapse Q = ln(ae _12<TV ), we find that the requirement for 
well-posedness D = A = 12cra > yields a necessary 
condition for a, which is of course the same as (|46l) . 

To conclude this sub-section we note that the results of 
the analysis presented above do not depend on the order 
of linearization and enforcement of algebraic constraints 
of BSSN. Instead of linearizing the unconstrained BSSN 
first, one could have chosen to reduce the number of vari- 
ables of BSSN, by eliminating as many variables as there 
are algebraic constraints and then linearize the reduced 
system. It is straightforward to check that both ways 
lead exactly to the same minimal sets. 
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C. BSSN and Differential Gauges 

One can perform the same constrained perturbation 
analysis for the BSSN formulation in conjunction with 
the differential gauges considered in the ADM analysis 
above. However, the numerical relativity community has 
recently been resorting to non-trivial shift conditions, for 
example [13, HH , as an attempt to accurately evolve black 
hole binaries. Therefore, instead of analyzing the same 
gauges as in the ADM analysis above, here we focus on 
the elliptic "Gamma freezing" condition [24|, which was 
formulated for the BSSN formulation 



&r = o 



djd t f j = o. 



(112) 



This condition obviously " freezes" the evolution of the T z 
variables, hence the name. By use of (|112[) and the evolu- 
tion equations of P , from (|8"2")l , one obtains the following 
elliptic equations which the shift vector has to satisfy. 



2{d 3 a)A^ + 2a(r jk A k i - (2/3)^' {djK) 
•6A«(0 iV O)-0i(/^(&7 <J ') 



(113) 



(2/3) 7 ij ' (d k p k )) = 0. 



Here we consider only a ID perturbation approach about 
flat space to show that this gauge is good at least in the 
case considered. 

To reduce the system to first order form we define the 
derivatives of the shift vector as new variables Bf = 
dj/3 l . The derivatives of the shift satisfy (pS")) as in the 
ADM analysis. Of course the introduction of those 9 new 
variables leads to the introduction of 9 new constraints 
which have to be fully imposed. In the high frequency 
limit of small amplitude perturbations those constraints 
read djSj3 l — together with ([57]) . In terms of the new 
variables the linearized principal part of (|113[) becomes 
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k + \l l3 d k 8B 3 k - 2a\fidj5K + pdjST* = 0. 

(114) 

Equation (| 1 14|) can be treated as 3 more constraints in 
our approach, allowing us hence to eliminate the three in- 
dependent perturbations 5B^ . Thus, after imposing all 
available constraints the minimal set consists of 8 equa- 
tions, as it was expected. To simplify the analysis fur- 
ther we study this shift condition in conjunction with a 
lapse function of the form a = a(ip), that is the lapse 
depends only on the determinant of the 3 metric. This 
simplification results in the evolution equations of SA^j 
in (|106p - (|108[) being the same, but where A u = 0, so we 
will not write them here. However, the evolution equa- 
tions of SDij.i in (|106[) - (|108p are different, since we have 
to deal with the shift terms now. Those equations change 
as it is dictated by the linearized evolution equations of 



8Dij;k in (I92p . which in new variables read 



d t SD 



(■••) + 



_ dSB^ 



2, dSB 1 
o7« 



(115) 



Where (...) stands for the RHS of the evolution equa- 
tions of 5D ij; i in (fIDo ]) -([in5] l . All other terms in 1(52]) 
contribute to the low-order part only. If one imposes 
all possible constraints then the evolution equations of 
5Dij-i yield 

d t Dxi-i =0 
dtDisa =0 

<9 t -D 2 2;l 

d t D 2 3-i =(-..) 



l fll 0£>n ; i , dA u 



(116) 



The Jordan decomposition of the resulting system shows 
that the Jordan matrix is diagonal with the following 
eigenvalues on the diagonal {0, 0, 0, ft 1 , (3 1 +a, /3 1 — a, + 
a,^ 1 — a}. All eigenvalues are real and hence the sys- 
tem is well behaved. Surprisingly, the eigenvalues which 
correspond to the functional form of the lapse, that is, 
±ojA/6 are missing. Therefore, the "D-freezing" condi- 
tion gives rise to a well behaved constrained ID evolution, 
even if the lapse chosen is fixed or calculated by the maxi- 
mal slicing condition, because in ID perturbations about 
flat space this shift condition eliminates the gauge waves 
which correspond to the algebraic lapse considered when 
the shift vector is fixed. This is what one truly obtains, 
if one carries out the analysis of the " T-freezing" shift in 
conjunction with a fixed lapse or maximal slicing. 

Those results constitute a good indication of the well- 
posedness of the constrained BSSN evolution with this 
gauge condition. However, a definite answer requires a 
complete analysis, that is, consideration of planar per- 
turbations about an arbitrary spacetime. This is a very 
complicated task. We will address this in a future pa- 
per along with the other popular shift conditions, the 
parabolic and hyperbolic T- Driver" conditions [24[ in 
conjunction with the Bona-Maso family of slicing condi- 
tions. 



D. Equivalence of conditions for well-posedness of 
constrained evolution of ADM and BSSN 

We will now demonstrate that, in the limit of high 
frequency perturbations about flat space, the ADM con- 
dition (13511 which can be written as 



<91na 
A = — ViVj > 



(117) 



is equivalent to that of BSSN, equation (|104[) . We remind 
that 



%j = e 4 %j and tp = i In | 7 | 
Now using that dj = r yj lJ d'jij, we obtain 
^ 1 u 



fyij 12 



(118) 



(119) 
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and 



Since, 



<91n a <91n a dip d\n a d^kr, 



d-fij dip djij djkm d~/ i: 



(120) 



(121) 



we obtain: 
Sin a 



A — -y ij + e~ 4v A lj - e - 4 v> A fem — 7 4i . (122) 
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However, we are considering flat space so the latter be- 
comes 



<91n a . 1 
= A— 

5 7 » 12 



= A^6 lj + A ,J ' - A fero — 5 ij . (123) 



Therefore, if we consider that the vector along which we 
perturb is unit, the strong hyperbolicity condition of the 
ADM formulation yields: 

' Ua -ViVi = — (A + 12A^v iVj - 4A fcm 4 m ) > (124) 
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This last one is the same as condition (11091) . 

If one considers gauges for which the lapse function 
depends on the determinant of the three-metric and/or 
spacetime coordinates, as we have already shown the con- 
dition for well-posedness, reduces to 

din a 1 . 1 d\na „ 

A=- = — A = — > 0. (125) 

<91n 7 12 12 dip v ' 

which for both the ADM and BSSN language depends on 
the functional form of the lapse only. 



V. DISCUSSION AND CONCLUSIONS 

We have presented a general theory to study the well- 
posedness of constrained evolution of systems of quasi 
linear partial differential equations, which consist of a 
set of n first order evolution equations and a set of m 
first order constraint equations with m < n. We applied 
this theory to constrained evolution of 3+1 formulations 
of GR. In our analysis we explicitly took into account 
the Hamiltonian and momentum constraints as well as 
possible algebraic constraints on the evolution of high- 
frequency perturbations of solutions of Einstein's equa- 
tions. Our analysis revealed the existence of subsets of 
the linearized Einstein's equations that control the well- 
posedness of constrained evolution. 

We demonstrated that the well-posedness of ADM and 
3+1 formulations derived from ADM by adding combina- 
tions of constraints to the right-hand-side (RHS) of ADM 
and/or by a linear transformation of the dynamical ADM 
variables, depends entirely on the properties of the gauge 



and are equivalent to ADM on the surface on constraints. 
We note that our method concerns the constraint satis- 
fying modes only. Those are present in free evolution 
schemes, too. Therefore, a bad choice of gauge, which 
we define as one that produces ill-posed constrained evo- 
lutions, is bad for a free evolution, as well. However, a 
good choice of gauge for a constrained evolution scheme 
cannot in principle guarantee the well-posedness of a free 
evolution with the same gauge, due to the existence of 
constraint violating modes. 

Even on the surface of constraints we do not expect 
that all 3+1 formulation of GR which are derived from 
ADM by non-linear transformations and addition of ex- 
tra variables to have equivalent well-posedness properties 
when using the same gauges. For example the analy- 
sis of the exponential stretch rotation (ESR) formulation 
[26| , which is derived by a general non-linear exponential 
transformation of the ADM variables, shows that in the 
simplest case of geodesic slicing its behavior is elliptic, 
whereas ADM with the same gauge is weakly hyperbolic. 
The analysis of ESR will be the subject of a future paper. 

In this paper we also analyzed the BSSN 3+1 formu- 
lation which is derived by a non-linear transformation 
of the ADM variables and addition of extra dynamical 
variables. We were able to show that the well-posedness 
properties of BSSN and ADM on the surface of con- 
straints are similar for fixed and algebraic gauges. The 
results seem to indicate that, in general, the non-linear 
transformation of variables leading from ADM to BSSN 
does not change the well-posedness properties of the con- 
strained evolution when the same gauge is used. How- 
ever, the proof that the fully constrained evolution of 
BSSN is well-posed if and only if the constrained evolu- 
tion of ADM is well-posed, if such proof exists, is out of 
the scope of this paper. 

Our study shows that fixed gauges, that is, when 
the lapse function and the shift vector depend only on 
the spacetime coordinates, result in an ill-posed Cauchy 
problem for the constrained evolution of both ADM and 
BSSN as well as many other 3+1 formulations of GR. 
Algebraic gauges on the other hand can give rise to a 
well-posed constrained evolution provided that they sat- 
isfy (35]), (02]) or (fT04|) . In particular, fixed shift with 
the "harmonic" and "1 + log" slicing conditions, as well 
as with a densitized lapse having a > are all well 
behaved gauges. Our study of the Bona-Masso family 
of hyperbolic slicing conditions showed that it provides 
us with a well-posed constrained evolution. The study 
of well-posedness of constrained evolution with maximal 
slicing and fixed shift shows that it depends on the way 
this gauge is implemented. The algebraic implementa- 
tion j^Kij = leads to a well-posed evolution whereas 
the often used differential implementation (|4"5|) or ([50]) 
is ill-posed. The parabolic extension of maximal slicing 
with fixed shift leads to an ill-posed evolution. Finally, 
we demonstrated evidence that the constrained evolution 
in conjunction with the 'T- freezing" shift condition and 
an algebraic lapse leads to a well behaved constrained 
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evolution at least in the case of ID perturbations about 
flat space. However, a complete well-posedness analy- 
sis is still required and this will be a subject of a future 
paper. 

Our analysis demonstrates that the weak hyperbolicity 
associated with fixed gauges is directly related to the in- 
ertial deformation of a synchronous reference frame with 
time which, in a general non-linear case when the pertur- 
bations are not small, leads to the formation of caustics 
(see equation (|42j) and discussion following it). 

Finally, we note that gauge stability may be investi- 
gated more directly by considering variations of gauge 
degrees of freedom only. In general, this requires the 
analysis of a system of eight quasi-linear partial differ- 
ential equations presented in [l3j |. The main advan- 
tage of the method outlined in this paper is that, in 
addition to gauge conditions, it provides us with sub- 
sets which control the constrained evolution of spacetime. 
The method is also able to provide sufficient conditions of 
well-posedness, whereas the analysis in [l3[ gives only the 
necessary conditions. The subsets controlling the con- 
strained evolution can be used for construction of stable 
numerical schemes for 3+1 formulations of GR. This will 
be the subject of our future paper. 
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APPENDIX A: STRONG HYPERBOLICITY FOR 
CONSTRAINED EVOLUTION 

In this appendix we show that for uniformly diagonaliz- 
able systems there exists a uniformly bounded similarity 
transformation S which diagonalizes the principal matrix 
A, provided that the following two conditions are met: 1) 
A has real eigenvalues, which are analytic functions of 
and 2) The elements of A can be represented as ratios of 
analytic functions of Vk ■ 

To study well-posedness of a constrained evolution as 
outlined in this paper we use the m linearized constraint 
equations (|L3[) to eliminate some of the dynamical vari- 
ables. Equations (|13p are an under-determined set of m 
algebraic equations for the spatial derivatives of the n 
unknown variables, which in general can be solved for m 
of the n spatial derivatives of variables u. We can write 



equations ([13]) schematically as follows 
. du m 



dX 



( du i\ n 

K^) = o, 



(Al) 



where column vector of to of the n dynamical 

variables of the formulation, matrix C is m x to and de- 
pends on the direction Vi along which we perturb, and 
F m (vi, is a column vector with m components which 
are functions of the direction Vi and the spatial deriva- 
tives of the q = n — m dynamical variables left, and solve 
them as 



du r . 



dX 



-(C)- l F m (i 



dUq 



)■ 



(A2) 



Substitution of (|A2|) in equations (|12p leads to a set of 
q = n — m linear partial differential equations for q of the 
initial n variables, which is schematically given by l|14p. 
This elimination process includes inversion of C(vi), and 
thus division by its determinant, which is a polynomial 
in the components of the unit one form Uj. This may 
not be possible for every direction because there may be 
directions which make matrix C singular. In order to 
obtain the minimal set, the determinant \C\ has to be 
non- vanishing. The domain of a minimal set consists of 
all directions for which \C\ ^ 0. For singular directions 
v s we must use another set of to dynamical variables for 
which |C| 7^ in the singular direction v s . It can be 
shown that this is always possible for the GR equations. 

A transformation matrix S(vi) which diagonalizes A 
in (|14[) has the same domain as A and has non-zero de- 
terminant in its domain. However, when we approach 
a singular direction, the determinant of S or its inverse 
may tend to zero (or infinity) and then ^ may not be 
satisfied. However, the choice of eigenvectors and the 
corresponding transformation matrix S are not unique. 
The eigenvectors can be rescaled and this will change S. 

The systems analyzed in this paper for which the eigen- 
values are real and for which there exists a complete set 
of eigenvectors for all directions, we find that the eigen- 
values are analytic functions of (see (|33|) . (|34|) . ([76])'). 
For such systems we show below that it is always possible 
to rescale the eigenvectors in such a way that all rescaled 
eigenvectors will be analytic functions of Vf.- Then, ac- 
cording to 23] the transformation matrix S will satisfy 
and thus the system will be strongly hyperbolic and 
by definition well-posed. 

First consider matrix A. Its coefficients may be ratios 
of polynomials due to substitution of constraints. This 
is the case with the Einstein equations and all gauges we 
have studied in this paper. We write this schematically 
as 



(Mvk)) i 



ij qij(Vk) 



Pll 


Pig 


9ll 


qi q 


Pql 


Pll 


9,1 


Qqq 



(A3) 



where Pij(vk) and qij(vk) are polynomial (and hence an- 
alytic) functions of Vk- We further assume that A q has 
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real eigenvalues and a complete set of eigenvectors for all 
possible directions Vk in its domain. This is the case for 
algebraic gauges and the Bona-Masso hyperbolic gauges 
considered in this paper. 

Let Vi be a set of eigenvectors corresponding to the 
eigenvalues oji . They satisfy 



(A q - u^Vi = Aiivk)^ = 0, 



(A4) 



where I q is the gxg identity matrix. The coefficients of 
newly defined matrices Ai are ratios of analytic functions 
because the eigenvalues of A q are analytic functions. 



qi 3 {vk) 



Ell 






Qlq 






Iql 


1m 



(A5) 



Consider now a particular eigenvalue, e.g., u>% and a par- 
ticular eigenvector V\. Dropping the subscript 1 we can 
write (IA4I) as 



Ell y + my 

911 1 912 z 



Pi 2 



(A6) 



9,j2 



Iqq 



-V a = 



here subscripts of V indicate components of a particular 
eigenvector V = Vi we consider. If the rank of the ma- 
trix .A; is r, where r < q then there are only r linearly 
independent equations in (|A6[) . This means that we can 
solve the algebraic set for r of the q components of the 
eigenvector V, which we denote by V r . The remaining 
s = q — r components of V form a vector V 8 . We write 
the reduced system schematically as 



D r V r = D S V S = B r 



(A7) 



where D r is a r x r matrix and D s is a r x s matrix, 
which are both formed by elements of Ai- Therefore, 
both D r and D s have elements which are ratios of an- 
alytic functions of Vk- The non-zero components of V s 



can be chosen freely, for example as constants, so that 
components of B r will be ratios of analytic functions in 
Vk- By Cramer's rule the solution of system (|A7|) for the 
r unknown eigenvector components is 



(Yr)i 



I (AO. 
\DJ 



(A8) 



where {V r )i is the i-th component of V r and (D r )i is the 
matrix formed by replacing the i-th column of D r by 
the column vector B r . Therefore, all components of the 
eigenvector V can be expressed as ratios of analytic func- 
tions, since the sums and products of analytic functions 
are analytic, which we write schematically as 



Pi(Vk) 



(A9) 



Qi{v k ) 

If we multiply the eigenvector V by the product Q = 
TliQi(vk) the rescaled eigenvector will be an analytic 
function of Vk- We can carry out such rescaling for all 
eigenvectors. Then according to 23] (|A1[) is strongly hy- 
perbolic because A has real eigenvalues and a complete 
set of eigenvectors which are analytic functions of Vk for 
any direction. 

We now illustrate an example of rescaling that leads 
to a uniformly bounded S. Consider perturbations in 
the x 1 ^ 2 plane about flat spacetime for the ADM + 
densitized lapse system (a — C{xi) r y a ), for which the 
components V\ and Vi , of the unit one-form along which 
we perturb, are assumed non- vanishing. For this simple 
case q=30-22=8. The linearized evolution equations of a 
minimal set are in matrix notation: 



du du 

at = 8 9A' 



(A10) 



where 



U 



and 
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where u T is the transpose of u. When a > 0, this matrix always has real eigenvalues 
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{O, 0, — a, a, — a, a, — a\/2a, a\^2a} and complete set of 
eigenvectors in its domain, i.e. v% ^ 0, i>2 ^ 0. This 
is also the domain of the matrix (5) of eigenvectors of 
A$ as columns, which is a matrix that diagonalizes 
via a similarity transformation. The determinant of this 
matrix turns out to be: 



APPENDIX B: ILL-POSEDNESS OF PARABOLIC 
MAXIMAL SLICING 



\S\ 



S2a 



Vi 2 V2 5 



(A13) 



As we approach the singular points V\ = 0, v 2 = this 
determinant blows up and one cannot obtain an upper 
bound to satisfy (0). However, if we define a new matrix 
(within the same domain) by 



the determinant of this new matrix is 



\S\ 



'2a. 



(A14) 



(A15) 



This is a well behaved non-singular transformation that 
satisfies dU). 



In this appendix we demonstrate that the ADM con- 
strained evolution with parabolic maximal slicing is ill- 
posed. Since the parabolic maximal slicing is a second 
order differential equation we introduce new variables to 
achieve first order form at least to the evolution equa- 
tions of those variables, whose solution may determine 
the solution of the entire system. Let Ak = J^. The 
evolution equations for these variables can be easily ob- 
tained by commuting a time with a space derivative and 
after we perturb along a specified direction about a base 
solution and make use of the definitions of these vari- 
ables, the equations that describe the evolution of high 
frequency perturbations of Ak are (assuming without any 
loss of generality that the component of v along x 1 is 
non-zero): 



dSA k 
dt 



1 d 2 6A 1 
—Vk- 

Vl 



dX 2 



^ v k v£ dSAi 

ij ' 



Vl 



dX 



M 
dl 



8D 

dx 



dD 



dX 



dX 



d5K, 
dX 



(Bl) 



Now define JCij = Kij + eviViVjOt. Then equations (|59p ~ 
(j62[) yield for the evolution of high frequency perturba- 



tions of a,K,ij and D 



ij;k- 

dSa _ 1 dSAi 
~dt ~ 



evi 



dl 



dSICij 
dt 
dSDjj.k 
dt 



--aRh 



dStCij 
dl 



(B2) 



(B3) 



The momentum constraints (|23j) have exactly the same 
form for this newly defined "extrinsic curvature" variable 
Kij in the limit of high frequency perturbations, i.e. 



System (|Blj) - (|B3|) has now decoupled in two parts. One 
consists of the evolution equations for the A k variables 
and the lapse function, equations (|B1|) and (|B2[) . The 
other part consists of the evolution equations for Dy ; ft 



and K-ij, equations (IB3j) . governed by the exact same 
constraint equations as variables D^-k and Kij. The con- 
strained evolution of this last set has been analyzed in the 
section of fixed gauges and ADM and has been found to 
be ill-posed. Since the constrained solution of (|B3jl de- 
fines the solution of (|B1[) and (|B2[) . the entire system is 
ill-posed. And hence the constrained evolution with the 
parabolic extension of maximal slicing is weakly hyper- 
bolic and hence is ill-posed. 



l ms v. 



d5K„ 



, d5K r 



dX 



= 0. 



(B4) 
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